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Abstract. We present a new algorithm which is named the Dynamical 
Functional Particle Method, DFPM. It is based on the idea of formu- 
lating a finite dimensional damped dynamical system whose stationary 
points are the solution to the original equations. The resulting Hamil- 
tonian dynamical system makes it possible to apply efficient symplectic 
integrators. Other attractive properties of DFPM are that it has an 
exponential convergence rate, automatically includes a sparse formula- 
tion and in many cases can solve nonlinear problems without any special 
treatment. We study the convergence and convergence rate of DFPM. 
It is shown that for the discretized symmetric eigenvalue problems the 
computational complexity is given by O {n ( d+1 )/ d ^ , where d is the di- 
mension of the problem and N is the vector size. An illustrative example 
of this is made for the 2-dimensional Schrodinger equation. Comparisons 
are made with the standard numerical libraries ARPACK and LAPACK. 
The conjugated gradient method and shifted power method are tested 
as well. It is concluded that DFPM is both versatile and efficient. 



1. Introduction 



1.1. The Dynamical Functional Particle Method. The goal of this pa- 
per is to present an idea for solving equations by formulating a dynamical 
system whose stationary solution is the solution of the original equations. 
Examples of equations that can be solved are ordinary and partial differ- 
ential equations, linear or nonlinear system of equations, and particularly 
eigenvalue problems. In this section we begin by formulating the equation 
and the dynamical system in an abstract setting. We then give the corre- 
sponding finite dimensional formulation by discretizing the infinite dimen- 
sional problem. This discretized problem will then be analyzed and studied 
throughout the paper. 

Let J- be an operator and v = v(x), v : X — > M. k , fcsN, where X is a Banach 
space that will be defined by the actual problem setting. We consider the 
abstract equation 

(1.1) ?(v) = 

that could be, e.g., a differential equation. Further, a parameter t is in- 
troduced, interpreted as artificial time, which belongs to the interval T = 
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[to,oo). A related equation in u = u(x, t), u : X x T -> R k is formulated as 
(1.2) iJLUtt + r]Ut= F(u). 

The parameters /i = fj,(x,u(x,t),t), r] = r](x,u(x,t),t) are the mass and 
damping parameters. The idea in the infinite dimensional setting is to solve 
( |1.1[ ) by solving (1.2) in such a way that Ut, Utt — > when t —¥ t±, t\ < oo, 
i.e., limt^.t 1 u(x,t) = v(x). In addition, the two initial conditions u(to) and 
Ut(to) are applied. 



Both (1.1) and (1.2) need to be discretized to attain a numerical solution. 



For simplicity, we exemplify by applying finite differences but it is possible 
to use, e.g., finite elements, basis sets or any other method of discretization. 
We define a grid x%, X2, ■ ■ • and approximate v(xi) by v% and assume that the 



discretized version of (1.1) can be written as 
(1.3) 

where F~ : R n -> E 



Fi(vi ...,v n ) = 0,i 



1, 



Turning to the dynamical system (1.2) it is discretized such that Ui(t) ap- 
proximates u(xi,t) and (J,i(t) = fj,(xi,Ui(i),t), rji(t) = rj(xi,Ui(t),t) for i = 
1, ... ,n. Further, J~(u) is discretized as J~{v) in (1.3) and we approximate 



,n. 



(|L2| with the system of ordinary differential equations 
(1.4) HiVn + ruiii = Fi(uu ■ ■ -,u n ), i = 1, ... 

with initial conditions Ui(to), Ui(to). Our idea in the discrete setting is to 
solve (1.3) by solving (1.4) such that iii(t), Ui(t) — > when t — > t\, t\ < oo, 
i.e., lim^^ Ui(t) = Vi. The overall approach for solving (1.1) using (1.4) is 
named the Dynamical Functional Particle Method, DFPM. 



1.2. Related work and topics. In a recent mechanics oriented article the 
connection between classical particle methods and differential equations were 
studied |19J . This work had a clear focus on the physical understanding 
and mechanical properties. The present work, however, turns the focus to- 
wards the mathematical aspects in the attempt to answer questions related 
to convergence, rate of convergence and the underlying reasons why DFPM 
is seen to be efficient for some mathematical problems. The idea of study- 
ing dynamical particle systems certainly has its origin in basic physics and 
astronomy. The assumption there is that all matter consists of particles. 
Their interactions are known and they follow the equations of motion. The 
basic idea DFPM, however, is that the "forces" and "particles" instead are 
viewed as mathematical abstract objects rather than physical. Thus mathe- 
matical quasi-particles are formed. Their interactions are determined by the 
functional equation at hand. From the mechanical point of view the quasi 
particles in (1.3) have masses jXi and all follow a dissipated motion governed 
by r/j. Such Hamiltonian systems have many interesting and useful proper- 
ties, see, e.g., |21j . In Hamiltonian systems it is well known that symplectic 
integration techniques are especially attractive |30[ [25] . 



The idea of solving a time dependent problem to get the stationary solution 
has also previously been applied in mathematics. A simple example is the 
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solution of an elliptic PDE such as in the heat equation, see Sincovec and 
Madsen |32j . Indeed, steady state solutions are often the objective when sim- 
ulating time-dependent PDE systems. Since the stationary state is seeked, 
the evolution of the system is considered to take place in artificial time. The 
concept of artificial time is further discussed and analyzed in [6]. 

A general approach is that of continuation, see [3] for an introduction, where 
is embedded in a family of problems depending on a parameter s, i.e., 

(1.5) T(u;s) = 



where J-(u; 0) = F{u) = 0. Thus, solving (1.5) for s = is equivalent to 



solving (1.1) and it is assumed that solving (1.5) for some s, say s = 1, is 



computationally cheap. The solution to (1.5) is found by solving a sequence 
of problems for values of s decreasing from 1 to 0. A general package for 
continuation methods with additional references may be found in Watson 
et al. |3l]. Further, see Nocedal and Wright |27] for a discussion in the 
context of optimization and Ascher, Mattheij and Russell [7j for boundary 
value ODEs. DFPM can in principle be viewed as a sub-method to the 
group of continuation methods. However, as far as the authors know, the 
concrete application of a second order system (Hamiltonian dynamics) to 
solve equations and the corresponding analysis as presented here is novel. 



Other works where (1.2 ) appear are for example the damped harmonic oscil- 
lator in classical mechanics, the damped wave equation, |28| and the heavy 
ball with friction j5] . These problem settings are specific examples of physical 
systems and not developed to solve equations in general. 

In |13| [T3] iterative processes to solve, e.g., eigenvalue problems are con- 
sidered as (gradient driven) dynamical systems. So called fictitious time is 
used in |33] where, e.g., Dirichlet boundary value problem of quasilinear el- 
liptic equation is numerically solved by using the concept of fictitious time 
integration method. The inverse problem of recovering a distributed param- 
eter model is considered in [6] using the first order ODE attained from the 
necessary optimality conditions of the inverse problem. 

First order systems, mainly in the form ut = J~(u), have been used to solve 
different kinds of equations J-(v ) = 0, both as a general approach and in- 
tended for specific mathematical problems. It is of interest to briefly con- 
sider the difference between the first order differential equation Ut = J~{u) 
and the second order approach, DFPM. Suppose for simplicity that a dis- 
cretization is made by finite differences leading to a system of equations 
u% = Fi(ui, . . . ,u n , Xi). Consider an example where the functional J-(u) 
contains a derivative w.r.t. u (A) and other functions of u and x (B). Then 
Fi(u\, . . . ,u n , Xi) = Ai(u\, . . . ,u n , Xi) / h p + Bi(iti, Xi) . Dimensional analysis 
then gives that 

M = M = [Mm, ...,u n , Xi )/hP] = M = M 

Given a certain component Fi we have that Xi is not variable, so [A] = 
[u]. We see that the dimension of time is related to the discretization, i.e., 
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t = O (h p ). In a similar way it can be shown that for the second order 
differential equation (DFPM) one instead have that t = O (/i P//2 ). In a 
numerical integration, the dimension of t must still be valid. This means 
that also the maximum timestep (for which a numerical algorithm is stable) 
is given by At max = O (h p ) for the first order case, and Atmax = O (hP/ 2 ) 
for the second order equation. Consider for example the case of a central 
difference formula, i.e., p=2. For the first order differential equation we see 
that At m ax will have to be very small as finer meshes are selected. This will 
lead to a less efficient computational complexity for the first order system. 
The complexity of DFPM will be further discussed in section [HJ 

As we shall see, the DFPM algorithm seems attractive due to several rea- 
sons. The most interesting points that will be studied are related to compu- 
tational complexity, easiness of implementation, Hamiltonian dynamics and 
its relation to the total evolution time, stability and cheapness of symplectic 
integration, exponential convergence and the existence of potential energy. 



1.3. The outline of the paper. The outline of the paper is based on il- 
lustrating the versatility of DFPM and to analyze the convergence aspects 
of the dynamical system (1.4). In order to introduce the reader to DFPM, 
the damped harmonic oscillator is revisited. DFPM clearly has a close re- 
lationship to this type of classical system. It is then important to remind 
the reader of the close connection to Hamiltonian dynamics in particular 
the existence of a potential function. In such a case where the functional 
is conservative any extreme value of its potential function is a solution to 
the original problem (1.3). Specifically if the potential has a minimum the 
solution of (1.4) will converge asymptotically. This is dealt with in Section [4] 
A Lyapunov function is applied to show asymptotic convergence in Section 



4.1 In Section 4.2 we analyze the linearization of DFPM and give precise 



statements for local asymptotic convergence valid close to the solution and 
for linear problems such as systems of linear equations. The rate of con- 
vergence is treated in Section [5] with four subsections treating the general 
problem (1.4) when there is a Lyapunov function, the linearized problem, the 
choice of damping, and examples, respectively. In Section 5.1 the Lyapunov 
function is used to state a general theorem that together with additional as- 
sumptions on the Lyapunov function gives an exponential convergence rate. 
This theorem is then specialized to the case when there exists a potential. 
The linearized problem is analyzed in Section 5.2 where we first treat the 
case with one scalar damping and then discuss the possibility to choose an 
optimal general damping matrix. The conclusions drawn from the choice 
of damping in the linear case are used in Section 5.3 to formulate a local 
strategy for the choice of optimal damping. To demonstrate the efficiency 
of DFPM we report in the end of the article several examples. The most 
noteworthy is the efficiency for treating symmetric eigenvalue problems. It is 
shown that DFPM is order of magnitudes faster than the standard software 
ARPACK p]. Finally, in Section we make some conclusions, discuss open 
problems as well as suggestions for future works. 
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2. TWO ILLUMINATING EXAMPLES 



2.1. The damped harmonic oscillator. Despite its triviality an impor- 
tant example related to DFPM is the damped harmonic oscillator where 
F(u) = — ku, k > which we include for later reference since it illustrates 



many properties of (1.4). In this case the equation at hand is given by 



(2.1) {iu + r]u = —ku. 

In DFPM, as well as here, we set the inital condition ii (0) = 0. The initial 
condition for u may be set arbitrary. In mechanics the parameters \x > 0, 
r] > and k > correspond to particle mass, damping constant and spring 
constant, respectively. The time-dependent solution is given explicitly by 



u = cie" 71 * + c 2 e" 72 * 



where 

(2.2) 71,2 




k 

and ci, C2 are constants given by the initial conditions. Although in mechan- 



ics it is clear that all parameters in (2.1) are physically positive, it is worth 



while to make a comment why this is so. Consider the case /J, < 0, f] < 0. 
The roots are then real with one positive and one negative root so u (t) will 
diverge. The situation for fi < 0, i] > is similar with one positive real root 
and no convergence. When fi > 0, rj < the roots may be complex but one 
root will always have a positive real part and the solution will not converge. 
Thus, the only possible choice for convergence into a stationary solution is 
to apply positive parameters. 



There are three different regimes of the parameters that will effect the con- 
vergence: the under critical damping, r] < 2^/kJl which shows an oscillatory 
convergence, the critical damping, r/ = 2\/kji giving exponential convergence, 
and over critical damping, ij > 2^fk~ii resulting in a slower exponential con- 
vergence. The critical damped system is known to be the fastest way for 
the system to return to its equilibrium (i.e., the stationary solution) [2]. It 
will be illustrative to return to this example later when considering various 
aspects of convergence and convergence rate of DFPM in Sections [4] and [5} 



2.2. A symmetric eigenvalue problem. Symmetric eigenvalue problems 
are of great importance in physics and technology. It is also a relatively 
straight forward example to illustrate how the DFPM algorithm is applied. 
Consider the eigenvalue problem 

(2.3) Aw = Av 

where A is a symmetric matrix with normalization ||v|| = 1. The DFPM 



equation (1.4) is not directly applicable since the eigenvalue A is unknown. 
However, it is well known, see [22], that the stationary solutions of the 
Rayleigh quotient 

p(v) = v T Av, ||v|| = 1 
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are eigenvectors of (2.3). Specifically, we have that 

argminp(v) = A 

T. 



where A m j n is the smallest eigenvalue of A. Thus, one way to formulate the 
functional vector F is 



F = -Au + (u T 4u) u, 1 1 u 1 1 = 1 
where u = u (t). The DFPM equation ( |1.4[ ) is then given by 

(2.4) Mu + Nil = -Au+ (u T ^u) u, ||u|| = 1 

where M = diag(/zi, . . . , /i n ), N = diag(ryi, . . . , rj n ). This procedure will 
yield Xmin and its corresponding eigenvector u. We shall see later that by 
replacing F with — F we instead get A max and its corresponding eigenvector. 
There are various strategies to get the other solutions. One possibility is 
to apply a Gram-Schmidt process [llj. Often in applications, only a few of 
the lowest solutions are of interest. The reader should note that in practice 
the matrix A never needs to be formulated explicitly. In DFPM one instead 
works with the components of the functional vector Fj. The mechanical 
interpretation is of course that this is the force acting on particle i. The 
formulation therefore automatically becomes sparse. This is later illustrated 
in Section |6l 



3. The Potential and Total Energy 



DFPM can be considered as a many-particle system in classical mechanics 
where F{ is the force acting on particle i, //j is its point mass, and —rjiiii is the 
damping force |21j . In this section we revisit the concept of a conservative 
force field and thus the existence of a many-particle potential. In DFPM the 
functional is not necessarily conservative, but if it is, the analysis is greatly 
simplified. By using the results in this section, we shall see in Section |4.1 



that for a convex potential, the stationary solution to (1.4) corresponds to a 
minimum of the many-particle potential. 



We start by taking the view that (1.3) is a vector field in 



p n . 



(3.1) F = (F 1 ,F 2 ,...,F n ), F i = F i { Ul ,u 2 ,...,u n ). 



Definition 1. The vector field F in (3.1) is conservative if there exists a 



potential function V : R n -> R such that F = - W. 



For any conservative field F we have that 

F(v) = 0« W(v) = 0. 
Thus, any solution to (1.3) is an extreme value of the potential V, i.e., 



a minimum, maximum or saddle point. In other words, 



solving (1.3) by 
We 



DFPM is equivalent to finding the extreme points of the potential V . 
will explore this fact further when analyzing the convergence of DFPM in 
Section [4] 
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By differentiating W and assuming that V is at least twice continuously 
differentiable, we get 



/ x dFi dF j 

3.2 W ± ~^ A = 0,l<i,j<n. 

OUj OUi 

as a necessary and sufficient condition for F to be conservative. Note that 



one good example fulfilling the condition (3.2) is the force field F(u) = Au 
where A is a symmetric matrix. We shall see later that this fact is very 
useful for symmetric eigenproblems. 



It is possible to derive the condition (3.2) that is interesting in its own since 



it contains the possibility to consider equations on manifolds. Consider the 
(work) 1-form ip = ^Fjduj, see |24j for a definition of fc-forms. The 1-form 
p is said to be closed if dip = and Poincares Lemma [15] implies that 
any closed form is exact, i.e., in our context has a potential, say V, that is 
dV = p. There is no ambiguity to say that V is a potential as in Definition 
[T] We have the following results for the vector field in (3.1). 



Theorem 2. The following statements are equivalent: 

1. F is conservative, that is F = — W 

2. dV = ip 

3. jy. F ■ dr is independent of the path V 
4- §p F • dr = for all closed paths T. 



Proof The equivalence of 1 and 2 is trivial. The equivalence of 1 and 3 can 
be found in [24] , Now, assume that 3. is true. Let p and q be two arbitrary 
points, and let 71 and 72 be two piecewise smooth paths from p to q. Define 
r as the closed path which first goes from p to q, via 71, and then from q to 
p via —72- Then, since f F • dr = f F-dr, we get that 



F • dr- / F • dr = / F • dr. 

71 •> 72 Jr 

Since p and q are arbitrary points, and 71 and 72 are arbitrary, the closed 
path r is arbitrary, and therefore the implication 3 to 4 is proved. The 
implication 4 to 3 is similar. □ 



If a potential exists it can be derived from the discretized equations by 
calculating the work, W, simply integrating along any path, say, from to 
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u. For example it is possible to use coordinate directions as 

Ul U2 

(3.3) W = j Fi(si,0,...,0)dsi + J F 2 (u 1 S2,0,...,0)ds 2 + ... 







+ J F n (u 1 ,u 2 , ■ ■ ■ ,u n -i,s n )ds n = -V. 



3.0.1. A revisit to the symmetric eigenvalue problem. Recall that DFPM 



equation for the symmetric eigenvalue problem (2.4). The corresponding 
vector field 

F = -Au + (u r Au) u, 1 1 u 1 1 = 1 

is conservative with the potential 

(3.4) V(u) = \^ Ta ^ ll u ll = 1 

To prove this it would at a first glance seem natural to find the gradi- 
ent of V(u). However, the normalization ||u|| = 1 complicates this some- 
what. This can be treated by investigate the gradient on the sphere S n ~ l = 
{u G W 1 :, ||u|| = 1}. Denote the tangent space to the sphere at a point u 
as T u (5 n_1 ). By using the Euclidean metric (the 2-norm), the gradient of 
V at u is the unique vector V S n~iV(u) G T u (S' n_1 ) such that 

W(u) T t = V S n-lV(\l) T t 

for all tangent vectors t G T u (S' Tl_1 ) where VV(u) is the usual gradient in 
M n . Solving this equation for Vgn-i V(u) by realizing that V^n-iV^u) is the 
projection of W on T u (5 n_1 ) we get 

V s ,-iF(u) = (J - n(u)n(u) T )W(u) = VV(u) - (n(u) T Vy(u)) n(u) 

where n(u) is the normal to T u (5' n ~ 1 ). Since, for S"™ -1 we have n(u) = u 
and we get 

V S n-iV(u) = Au - (u T ^u) u = -F(u) 
showing that V in fact is a potential to the vector field F. 



4. Asymptotic convergence 



In this section we investigate the convergence properties of the solution u(i) 
given by (1.4). Since we are interested in the asymptotic solution we will use 
stability theory for dynamical systems, see [26J, namely the use of a Lyapunov 
function and local linear analysis. However, it is generally difficult to find a 
Lyapunov function. Consequently, we start with the case where there exists 
a potential and where the Lyapunov function can be chosen as the total 
energy. If no potential exists we are left with the linear stability analysis in 
Section [O 



For simplicity and without loss of generality we may assume that ^ = 1 and 
that the solution of (1.3) is u = 0. 
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4.1. Using the Potential. In this section we assume that there exists a 



potential, V(u), to the given equations in (1.3). Then (1.4) may be written 
as 

(4.1) ii + iVu = -W 

where N = diag(r/i, . . . , r] n ). The energy functional (the Lyapunov function) 
is given by 

(4.2) E = T + V 
where 



It, 



2 ' 



is the kinetic energy. We then have the following important result to be used 
in the analysis of the asymptotic convergence analysis. 

Lemma 3. Assume that rji > 0, i = l,...,n. For the given solution of 



the dynamical system (1.4) the energy functional defined in (4-2) is a non- 
increasing function. 

Proof. From the definition of E in ( |4.2[ ) and W = — F we have 
dE dT dV v-^ dV v-^ 

and since F{ = Uj + rjiUi we get 

(4.3) f = -E^B0- 

i 

□ 

Lemma [3] tells us that the energy is non-increasing which is not surprising 
from a mechanical point of view since the damping will decrease the total 
amount of energy and there are no additional sources of energy in the system. 

The next theorem is taken from |26j . The proof is omitted. 
Theorem 4. Consider the autonomous system 

(4.4) w = G(w) 

where G : O — > M n is a continuous function defined on a domain in 
M n containing the origin and G(0) = 0. Assume that the Lyapunov function 
L is non-negative with respect to {4-4^ f or all w £ Q and such that for some 
constant c £ R the set H c is a closed and bounded component of the set 
{w £ Q : L(w) < c}. Let M be the largest invariant set in the set 

„ dL , . 
w G Q, : — (w) = 
dt y ' 



Then every solution w(i) of the autonomous system (4-4) with w(to) G H c 
approaches the set M as t — > oo. 



Using Theorem [4] we are now able to show our main result in this section for 
the asymptotic convergence of (4.1). 
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u, to (1.3) and a potential 
then for 



Theorem 5. Assume that there exists a solution, 
V to F. If V is globally convex, i.e., convex in the whole of 
any initial starting point in 
convergent to u. If V is locally convex in a neighbourhood T of u then for 
any initial starting point in T the solution of (4-1) will be asymptotically 
convergent to u. 



the solution of (4-1) will be asymptotically 



Proof. Rewrite DFPM (4.1) as a system by letting v = u as 

r\ U V 

^ ' [ v J ~ [ -jvv - v u y 

We want to use Theorem [4] to prove asymptotic convergence of DFPM. From 



Lemma [3] we have that dE / dt 
and only if v = 0. Define 



V iVv < and therefore dE/dt = if 



w 



u 

v 



Let M = Z = {u : dE/dt = 0} = {w : v = 0} be the invariant set in 
Theorem [4] Then, by Theorem [4] again, the solution w to the system (4.5) 
approaches the set M as t — > 00. If w £ M then v = u = 0, so from (4.5) 
we have that v = —V U V 7^0ifu^0 = u and w can not remain in the set 
M if u 7^ u. We need to verify that u = as t — > 00, but this follows from 
the fact that E is non-increasing. Hence w — )• as t — ► 00. □ 



4.2. Linear Convergence Analysis. Without a Lyapunov function and 
with no existing potential one is left with a local linear stability analysis. 



Such an analysis is based on the linearization of the dynamical system (1.4) 
at a solution , u = 0, to ( |1.3| ). Define J(u) as the Jacobian of F. From the 
Taylor expansion F(u) = F(0) + J(0)u + 0(||u 2 ||) we define the linearized 



problem to (1.4) as 

(4.6) Mii + Nu = F(0) + J(0)u 

where 

M=diag(£i, . . . , Am), iV=diag(j?i, ...,fj m; 



and jli, fji are the values of /i, r] at u. Thus, the local convergence can be an- 



alyzed by analyzing the linear system (4.6) which for notational convenience 
is written 

(4.7) Mu + Nu + Au = b 

where M, N, A £ M nxn ,b G W n . In [17J sufficient conditions are given for 



(4.7) to have asymptotically stable solutions. In order to state these condi- 
tions we need some additional notation. Consider the homogeneous equation, 
i.e., 

(4.8) Mii + Nu + Au = 0. 



By inserting the eigensolution e^*Vj into (4.8) we get the equation 



(4.9) (4 2 M + ^ + A) Vi = 
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for the eigenvalue £j and eigenvector Vj of A. Let Re(Vj) , Im(vj) denote 
the real and imaginary part of Vj, respectively. Introduce the symmetric and 
anti-symmetric parts of M, N, A as Ms, Ma, N$, ■ ■■■ and define 

Si(M) = Re(v i ) T M 5 Re(v i )+Im(v i ) T M 5 Im(v i ), ai (M) = 2Re(v i ) T M A Im(v J 

with corresponding definitions for Si(N),ai(N),Si(A),ai(A). By applying 
the general Hurwitz criterion, see e.g. [20] , to ( |4.9[ ) we get the following 
theorem and corollary. For a detailed presentation of the proofs we refer to 

mi- 



Theorem 6. The solution to (4.7) will converge asymptotically if and only 
if 

Si(M) Si (N) + ai (M)ai{N) > 

and 

( Si (N) Si (A) + ai (N)ai(A)) ( Si (M) Si (N) + ai (M) ai (N))- 

- (si(M)ai(A) - a l {M)s l {A) f > 
Corollary 7. If M, N are positive definite and A has eigenvalues with pos- 



itive real parts then the solution to (4.7) will converge asymptotically. 



Let us now return to the question of local convergence for (1.4) and state 
the following theorem. 

Theorem 8. Define J(u) as the Jacobian o/F. Assume that there exists a 



solution, u, to (1.3). Further, assume that 

M=diaq{jix, Am), N=diag(fji, ...,fj m ) 
where jXi,f]i are the values of fx,r] at u. Then DFPM will converge asymp- 



totically for any initial starting point of (1.4) close enough to u if and only if 
the conditions in Theorem^ are fulfilled where M = M , N = N, A = — J(u). 
Further, if M,N are positive definite and J(u) has eigenvalues with negative 
real parts then DFPM will converge asymptotically for any initial starting 



point of (1-4) close enough to u. 



Proof. The first statement in the theorem follows directly from Theorem [6] 
and the second from Corollary [7] □ 



5. Convergence rate 

5.1. General results on convergence rate. Sharp estimates of the con- 
vergence rate for DFPM in a general case is difficult and not realistic. How- 
ever, we shall give some important special cases that is relevant for solving 
equations with DFPM, i.e., to achieve fast exponential convergence. We em- 



phasize that this is crucial to attain an efficient method for solving ( 1.3 ). We 



again assume without loss of generality that the solution of (1.3) is u = 0. 
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Definition 9. The solution, u, to (1.3) is locally exponentially stable and 



the solution u(t) to (1.4) has a local exponential convergence rate if there 



exists an a > and for every e > and every to > 0, there exists a 5 (e) > 
such that for all solutions of (1.4) ||u(t)|| < ee~ a (*~*°) for all t > to whenever 
ll u (^o)|| < The solution, u, to (1.3) is globally exponentially stable 



and the solution u(t) to (1.4) has a global exponential convergence rate if 
for /3 > there exists k(j3) > such that ||u(i)|| < A;(/3)e~ a ( t_ *°) whenever 
||u(t )|| < k(/3). 



We begin by stating a general theorem from |26| giving one possible formula- 
tion of the requirements for exponential convergence based on the existence 
of a Lyapunov function L(u, u). 

Theorem 10. Assume that there exist a Lyapunov function L = L(u, u) : 

u ^ 
u 



B(r) -> M , B(r) 
such that 



< r > and four positive constants c±, C2, C3, p 



u 
u 



<L(u, U)<C2 



u 
u 



and 



for all 



dL 
~dl 



< 



-C3 



u 
u 



u 
u 



E B(r). Then the solution, u, to 



1.3) is locally exponentially 



stable and the solution u(t) to (I.4) has a loca l exponential convergence rate. 
If B(r) is the whole M 2n the solution u to (1.3) is globally exponentially stable 



and the solution u(t) to \1.$ has a global exponential convergence rate. 

In the case that there exists a potential we can choose the Lyapunov func- 
tion as the total energy E in (4.2). We will state a theorem that proves 



exponential convergence in this case that is slightly different from Theorem 

QZU 

Theorem 11. Assume that there exists a potential V(u) which is locally 
convex in a neighbourhood of the solution, u, to \1.3ty and satisfies the bound 

(5.1) Cl ||uf<y(u)<c 2 ||uf 

where c\ y Ci and p are positive constants. Then u is locally exponentially 



stable and the solution u(i) to (1-4) has a local exponential convergence rate. 
If the potential is convex and satisfies the bound (5.1) in the whole W 1 the 
solution u(t) to (1.4) has a global exponential convergence rate. 



Proof. From (5.1 ) and the definition of kinetic energy we have that the total 



energy is bounded as 
(5.2) 



c l II U H P + ^mai ||u|| 2 <E<C2 ||u|| P + ^fl m ax [|u| 



Thus, from Lemma [3] we have 



dE 

dt 



< 



U 
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If V is convex we see from Theorem [5] that ||u|| > unless at the solution u 
and therefore there exists a constant 7 > such that 

dE 



and thus E < £ e" 7(i ~ io) . From (5.2) we see that 



d ||uf + -/w ||u|| 2 < E < Eoe-^-^ 



and this implies 

for some positive constant c depending only on the initial conditions. □ 



1 I, ^ -±~{(t-t ) 
u < ce p 



5.2. Convergence rate for linear problems. Consider again the damped 



harmonic oscillator (2.1). Obviously, the convergence rate is exponential for 
all different dampings. However, the fastest convergence rate is achieved for 
the case where the damping is chosen as critical damping which can be seen 



in (2.2): The negative real part below critical damping is r]/2 and therefore 
r] should be as large as possible. However, as soon as critical damping is 
exceeded, one of the roots will be real and the negative real part will increase 
for the larger real root. 



This property is inherited for more general linear problems defined by (4.7) 
which we now shall investigate further. We will assume that M and N are 
diagonal matrices with positive elements and then we can, without loss of 
generality, assume that M = /, b = and consider the system 

(5.3) u + Nix + Aa = 0, iV=diag(r/i, . . . , r/ n ) 

where A £ M nxn . This linear system can be restated as 



u 







I ' 




u 


V 




-A 


-N 




v 



and since this is a linear autonomous system of first order differential equa- 
tions, it is well known that any convergent solution will have exponential 
convergence rate. However, we are interested in solving the linear equations 
as fast as possible and then it is reasonable to try to answer the question: 
How fast is the exponential rate of convergence? This is a difficult question 
for a general A and N because of the following result from linear system 
theory [29] . 

Theorem 12. Any two square matrices that are diagonalizable have the same 
similarity transformation if and only if they commute. 



For the problem (5.3), Theorem 12 means that AN = NA is a necessary 
and sufficient condition for having a similarity transformation such that 
A = TAaT' 1 ^ = TAjvT" 1 where A A ,A N are diagonal matrices. In other 
words, the two matrices A and iV has to commute in order to decouple the 



system (5.3) into n one dimensional damped harmonic oscillators where the 
optimal damping is critical damping as shown earlier. Note that the special 
case N = rjl where all damping parameters are the same, trivially commutes 
with any matrix A. 
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FIGURE 5.1. Smallest eigenvalue of the Hessian \7 2 V(u) for 
the potential (5.5 ). 



5.3. A discussion of optimal damping. The problem of choosing the 



damping in (1.4) such that the asymptotic solution is attained, to some 
precision, in minimal time is difficult, see, e.g., [12j,[31j for some special 



cases. Moreover, from a practical point of view this is not very interesting 
since it requires a priori knowledge of the solution. A more interesting 
approach is to choose the damping according to a local measure of curvature 
which we will discuss briefly. Assume that the solution to (1.3) is u = and 
that there exists a potential V that is convex with V(u) = and a positive 
definite Hessian V 2 V(u). Consider the case with a single damping parameter 
V = v(t)- Then a Taylor expansion at u in (1.4) gives the approximate 
problem 



(5.4) 



u + Tj(t)u = - V 2 V(0)u 



From the linear case treated in Section 4.2 the optimal damping for (5.4) 
is rj ~ 2y/ A m j n (V 2 y(0)) where X m in(') denotes the smallest positive eigen- 
value. Now, consider any u(t),t > to then we conjecture that a good choice 
of damping is rj(t) = 2y / A m j n (V 2 y(u(t))). To illustrate the possibilities of 
this choice of damping we given an example with a potential 



(5.5) 



V 



■2*1 



that is globally convex with a minimum at U\ = u% = 0. In Figure 5.1 the 
eigenvalue distribution is shown for the smallest of the eigenvalues of V 2 V. 
Indeed, looking at the trajectories in Figure 5.2 for the choice r](t) = 1 (solid 
line) and rj(t) = 1.9y A m j n (V 2 V (u(t))) (curve indicated with '*') it is clearly 
seen how the choice of damping affects the convergence. In fact, the effect 
is rather striking and further analysis and tests of our conjecture is of great 
interest in order to improve DFPM. 
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FIGURE 5.2. Trajectories for damping rj(t) = 1 (solid line) 
and rj(t) = 1.9 v / A min (V 2 y(u(t))) (curve indicated with '*'). 



6. DFPM EXAMPLE FOR THE HELIUM ATOM 



A relevant numerical application is to study the s-limit case of the Helium 
ground state energy. This example is often used in atomic physics literature 
as a benchmark to study numerical accuracy and efficiency The equation at 
hand is the Schrodinger equation. Due to electronic correlation and conse- 
quently discontinuities ("Cato cusps") this is often considered to be a tough 
problem. Another complication is the many-particle character leading to ex- 
tremely high dimensionality. The Helium example here only slightly touches 
these problems because the full correlation term has been neglected, i.e., 
that term is replaced by 1/max (ri,r2) resulting in only a 2D problem. This 
example is nevertherless sufficient to demonstrate many of the properties 
of DFPM. More complex examples have already been tested and DFPM 
remains relevant. However these fall outside the scope of the present work. 



Accordingly, consider the Schrodinger equation for the s-limit case of Helium 
US: 



Hv (n,r 2 ) 



(6.1) 



+ - 



1<P_ 
2drj 
1 



i d 2 



2drj 



2 

n 



>'2 



- + 



v (n, r 2 ) = Ev (r\, r 2 ) 



max (r\,ri) 

The boundary conditions are given by u(n,0) = v (0, r 2 ) = u (i?, r 2 ) = 
v (n, R) = 0. The discretization can be made by using central finite differ- 
ences with equidistant mesh sizes h = 0.1/l.l fc (for both Ari and Ar 2 ) where 
k is an integer chosen to get different problem sizes. The discretized version 
of Hv (ru, r 2 j) of a certain particle pij at the position (ru, r 2 j) becomes 



Hv (ni,r 2 j) w Hui 



2 h 2 



2Uj 



2ui 



+ 



Uj 



r 2 j max(r li ,r 2 j) 
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TABLE 1. Efficiency of DFPM for the s-limit Helium 1 S' groundstate. 



k 


N 


At 


Eo 


DFPM (s) 


DACG (s) 


ARPACK (s) 


S-Power (s) 


4 


23871 


0.066 


-2.863893321606(6) 


0.6 


0.9 


7.7 


14 


6 


34980 


0.055 


-2.868655504822(7) 


1.1 


1.5 


17 


34 


8 


51360 


0.045 


-2.871926990227(9) 


1.9 


2.6 


31 


68 


10 


75466 


0.037 


-2.874170330715(4) 


3.4 


4.7 


62 


153 


12 


110215 


0.031 


-2.875706726414(1) 


5.8 


8.2 


127 


319 


14 


161596 


0.026 


-2.876758055924(6) 


10 


14.4 


265 


676 


16 


237016 


0.021 


-2.877477040659(9) 


18 


24.9 


583 


1438 


18 


346528 


0.017 


-2.877968543434(3) 


33 


38.6 


1188 


3142 


20 


508536 


0.014 


-2.878304445684(1) 


58 


78.3 


2560 


6707 


22 


744810 


0.011 


-2.878533964475(9) 


103 


141.3 






24 


1090026 


0.009 


-2.878690772322(2) 


182 


249.7 






oo 


oo 




-2.8790287673(2) 











Note that the matrix H is never explicitly needed so the formulation is 
automatically sparse. The interaction functional component, i.e., the "force" 
acting on the particle pij at the position (ru,r2j) is given by Fij =< u\H\u > 
Uij—Huij . This can be compared with the equation ( 2.4 ) derived earlier. The 

notation < u\H\u > is the trapezoidal approximation to f Q R v (Hv) dr\dr2- 
The required norm ||f || = 1 is in the present context given by < u\u >= 1. 
The DFPM equation 1.4 for particle is thus given by 



(6.2) < u\H\u > Uij — Huij = fiilij + rjUij, < u\u >= 1 

In this case we apply constant mass and damping parameters (/i = 1 and 
7] = 1.54). The boundary is set to R = 15 which is sufficient for accurate 



ground state results. The integration method used to solve the ODE (6.2) 
is the symplectic Euler method, see e.g. [23J. The related Stormer-Verlet 
method was tested as well but gave no performance advantage. The test 
results are tabulated in Table 6.1. 

The first column shows k which determines the discretization h as mentioned 
above. In the second column, the corresponding total number of particles, 
N, is listed. Only a triangular domain needs to be computed due to even 
symmetry of the solution (i.e., v (r\,r2) = v (r2, ri)). Then the third column 
contains the maximum timestep At used in the Symplectic Euler method 
(depends on h). In the fourth column the eigenvalues E$ to 13 significant 
figures are listed. These values can easily be extrapolated to continuum (i.e., 
N —> oo). This extrapolated value is listed in the last line. In the final three 
columns we list the total CPU times in seconds to complete the computations 
to the desired accuracy. DFPM and three other methods are compared. 

A single C-code was written where the only difference was whether a function 
call was made to DFPM, DACG, ARPACK or S-Power. The DACG method 
(Deflation Accelerated Conjugated Gradient) is an advanced method to find 
some of the lower eigenvalues of a symmetric positive definite matrix. The 
algorithm is described in refs. [10} 19], Unfortunately the DACG algorithm 
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cannot be used directly because there are both negative and positive eigen- 



values present in the eigenvalue problem (6.1). Consequently, a shift of the 
potential (diagonal elements) had to be done. The shift for best performance 
was identified to be 3.9. The shifted power method (i.e., 'S-Power') is a sim- 
ple method that converges to the dominant eigenvalue [3J . This shift was also 
adjusted to get the maximum performance possible (about 0.5 (E max — Eq)). 
In the case of ARPACK, it is available at [I]. This iterative numerical library 
is based on the Implicitly Restarted Lanczos Method (IRLM). All tests were 
performed on a Linux PC with 1 GB primary memory and the CPU was a 
AMD Sempron 3600+ 2.0GHz. The compilers used were gcc and g77 version 
3.4.4. Both the numerical library and C-code were compiled using the op- 
timization: '-0'. The CPU times were measured using the Linux command 
'time' three times (the best reading is listed in Table 6.1). 

As can be seen in Table 6.1, the performance of DFPM is quite good consid- 
ering that no real optimization of DFPM has been attempted yet. DACG 
performs in parity with DFPM, although it is some 40% less efficient for this 
example. Both ARPACK and S-Power are far behind. However, ARPACK 
is seen to be more than twice as efficient as the S-Power method. For the 
smaller problems in Table 6.1, DFPM is one order of magnitude faster than 
ARPACK. For the larger problems, the advantage to DFPM is seen to ap- 
proach two orders of magnitudes. Also LAPACK was put to test (routine DS- 
BEVX). This routine computes selected eigenvalues and, optionally, eigen- 
vectors of a real symmetric band matrix. The parameter 'jobz' was set to 
compute eigenvalues only and the range was set to compute only the lowest 
eigenvalue. Unfortunately, due to how DSBEVX handles matrix memory 
allocation, large sizes, as applied here, is not realistic to compute. Even the 
smallest size N in Table 6.1 requires ~ 1000 s to complete. 

There are several reasons for good efficiency of DFPM. Firstly, the symplectic 
integration method for the second order differential equation allows a rela- 
tively large timestep. High numerical accuracy is not necessary during the 
evolution towards the stationary state. Only stability and speed is desired. 
Secondly, the cost per iteration is small because the symplectic integration 
method is as fast as Euler's method. Due to the damped dynamical sys- 
tem, the convergence rate is exponential in time. This is also theoretically 
consistent with Section [5] 

Further, in Fig. 6.1 it is seen that the computational complexity of ARPACK 
is given by O (-^V 2 ) . The DFPM complexity, however, is found to be O (Nzj. 

The cost of one iteration is the same for all the methods. Because of the 
sparsity of H, it is 0{N). The number of iterations is what separates the 
various algorithms. For ARPACK and S-Power the number of iterations 
are both 0{N). In fact, in the case of S-Power it is straight forward to 
prove that the number of iterations for a c?-dimensional problem is given by 
O (iV 2/d ), i.e., 0{N) for d = 2. DFPM and DACG, however, both show 

complexity O f° r the number of iterations. It is interesting to discuss 

this further for the DFPM algorithm. 
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FIGURE 6.1. Computational complexity of ARPACK and DFPM. 



Consider equation (6.2) and let us assume that after only a few iterations, 



< u\H\u >~ Eq. Most of the iterations are then spent solving: 



E Uij - Huij = fiilij + rjiiij 

for each one of the particles pij. The symplectic integration method is applied 
to approximate mj (t) on its way to the stationary solution. By applying a 
linear stability analysis of the symplectic method for the problem at hand, 
one finds that the maximum step size is 



\[E\ — Eq + -Eat-i — Eq 

where Ej^-i is the dominant eigenvalue. For the present problem, using the 
central difference formula, it is easy to show that this eigenvalue depends 
on the mesh size h according to, E^^\ = O (l//i 2 ). Since the mesh size 
only constitutes a minor correction to the lowest eigenvalues Eq and E\, 
we have that At max = O (h). The equidistant mesh size h and the total 
number of particles TV in a <i-dimensional problem are related, i.e. h ~ 
N^ 1 ^ d 1 thus Atmax = (A^™ 1 ^) . The number of iterations is therefore 
given by t/At max = O (A 1 /^), where t is the total evolution time until the 
stationary solution is achieved. The time t does not depend on N. That is, 
it is assumed that the mesh is fine enough and that the symplectic Euler 
algorithm approximately follows the true solution curve (i.e., the continuum 
case). The total complexity for DFPM is then given byO(A( (i+1 )/ d ). In the 
present test case, ri=2, and the behavior in Fig. 6.1 is thus explained. 

The presented benchmark indicates that DFPM is some 40% more efficient 
than DACG. DFPM can be further optimized by 1. allowing a varying 
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timestep during iterations, 2. preconditioning of the functional F and 3. 
allowing variation in the damping parameter during the iterations. However, 
also DACG can be optimized by a preconditioning of the matrix H. In the 
benchmark tests we experience that DFPM is more robust than DACG. 
DACG is quite sensitive to the selected shift. If the shift is small it has 
tendencies to diverge (despite that all eigenvalues are positive). If the shift is 
larger the convergence rate starts to suffer. The selected shift is |i?o|+l which 
gave the best convergence rate. Since DACG requires that the matrix H is 
positive definite, it would seem that DFPM is a better choice for Schrodinger 
problems where there often is a mix of positive and negative eigenvalues. 
Initially, the lowest (negative) eigenvalue is of course unknown which is why 
it can be difficult to apply DACG since one cannot assume a priori to know 
an appropriate shift. Another advantage is that the DFPM algorithm (1.4) 
remains the same also for nonlinear problems. The idea presented here is 
therefore quite versatile. The basic idea of DFPM as a dynamical system may 
also be attractive from a user's perspective since the algorithm is physically 
intuitive and very pedagogical. 



7. Conclusions and Future Work 



We have presented a new versatile method to solve equations that we call 
the Dynamical Functional Particle method, DFPM. It is based on the idea 
of formulating the equations as a finite dimensional damped dynamical sys- 
tem where the stationary points are the solution to the equations. The 
attractive properties of the method is that it has an exponential convergence 
rate, includes a sparse formulation for linear systems of equations and linear 
eigenvalue problems, and does not require an explicit solver for nonlinear 
problems. 

There are still a lot of interesting questions to be addressed. This includes the 



details for solving the ODE (1.4) in the most stable and efficient way. Moti- 
vated by the numerical tests reported in Section [6] we believe that symplectic 
solvers, see |5J, will be of great importance here. However, the stability prop- 
erties of the ODE solver is linked to the choice of parameters and especially 
the damping parameter. The key question is how to find the maximum time 
step retaining stability and getting a fast exponential convergence rate. We 
are currently working on these issues using the ideas presented here. 

DFPM has proved useful for very large sparse linear eigenvalue problems as 
indicated in Section [6] It is plausible that DFPM will also be efficient for 
very large and sparse linear system of equations. 

Since DFPM has a local exponential convergence rate it may be that it can 
be an alternative to the standard methods for nonlinear system of equations 
such as quasi-Newton or trust-region methods, see |16j . Moreover, if there 
exists a potential (or a Lyapunov function) DFPM has global convergence 
properties that can be useful for developing a solver for nonlinear problems 
with multiple solutions. 
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